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We have studied the static and dynamic magnetic properties of two-dimensional (2D) and quasi- 
two-dimensional, spin-S, quantum Heisenberg antiferromagnets (QHAF) diluted with spinless va- 
cancies. Using spin- wave theory and T-matrix approximation we have calculated the staggered 
magnetization, M(x,T), the neutron scattering dynamical structure factor, iS(k,u>), the 2D mag- 
netic correlation length, £(x,T), and the Neel temperature, Tn{x), for the quasi-2D case. We find 
that in 2D the hydrodynamic description of excitations in terms of spin-waves breaks down at the 
wavelength larger than tja ~ e*' Ax , x being impurity concentration and a the lattice spacing. We 
find the signatures of localization associated with the scale I and interpret it as the localization 
length of magnons. The spectral function for momenta a -1 3> k 3> t -1 consists of two distinct 
parts: (i) a damped quasiparticle peak at the energy cnk > u) S> ujq with abnormal damping 



Tfc ~ xcok, where u>o 



col. 



Co is the bare spin- wave velocity; (it) a non-Lorentian localization 



peak at u) ~ ton- For k < I' 1 these two structures merge and the spectrum becomes incoherent. 
The density of states acquires a constant term and exhibits an anomalous peak at w ~ uin associated 
with the low-energy localized excitations. These anomalies lead to a substantial enhancement of 
the magnetic specific heat, Cm, at low temperatures. Although the dynamical properties are signif- 
icantly modified we show that 2D is not the lower critical dimension for this problem. We find that 
at small x the average staggered magnetization at the magnetic site is M(x, 0) ~ S — A — B x, where 
A is the zero-point spin deviation and B ~ 0.21 is independent of the value of S; Neel temperature 
Tjv(x) ~ (1 - A s x) Tjv(O), where A a = n - 2/n + B/(S - A) is weakly S-dependent. Our results 
are in quantitative agreement with the recent Monte Carlo simulations and experimental data for 
S — 1/2, S = 1, and S = 5/2. In our approach long-range order persists up to a high concentration 
of impurities x c which is above the classical percolation threshold, x p « 0.41. This result suggests 
that long-range order is stable at small x and can be lost only around x ~ x p where approximations 
of our approach become invalid. 

PACS numbers: 71.10.-b, 75.10.Jm, 75.10.Nr, 75.40. Gb 



I. INTRODUCTION 



The discovery of superconducting cuprates (HTC) 
has motivated an enormous amount of studies in low- 
dimensional magnetic systems during the last fifteen 
years JO-p). Yet, the superconducting materials now 
form only a subfield in the activity of low-D quan- 
tum magnetism (see Ref. [0). The hope for an in- 
sight into the physics of HTC from the study of mag- 
netically related system has attracted much attention 
to the properties of diluted 2D, QHAF §-Q. One of 
such systems, La2Cui_ 2; Zn(Mg) a .04 (LCO), a quasi-2D, 
S = 1/2, QHAF diluted with spinless vacancies, has 
been a subject of great interest because of the possibil- 
ity of new quantum critical points (QCP) in its phase 
diagram. Earlier experimental data |17|] , while demon- 
strating that LCO shows much stronger stability against 
doping in comparison with the mobile hole doped com- 
pound La2- a: Sr a ;Cu04, have indicated the existence of a 
QCP at x = x c m 0.2, well below the classical perco- 
lation threshold. This finding was in a sharp contrast 
with the classical magnetic systems where dilution leads 
to the breaking of the magnetic bonds and long-range or- 
der (LRO) is lost only at the percolation threshold x p , a 



characteristic value of the dilution fraction x at which the 
last infinite cluster of connected spins disappears. For a 
2D square lattice x p rs 0.41 p8[. The existence of such 
a QCP below the percolation threshold was thought to 
be possible given the large amount of quantum fluctu- 
ations in the ground state of S = 1/2 system. How- 
ever, only recently the systematic experimental analysis 
of the diluted 2D AF has been performed in a wide range 
of doping p8|-p0| . Although there are several other ex- 
perimental realizations of a 2D QHAF on square lattice 
with S = 1/2 JH-H and S = 5/2 JUlf, the Cu °- 
based compounds are among the few which allow a di- 
rect probe via elastic and inelastic neutron scattering. At 
the same time, quantum Monte Carlo (MC) studies have 
provided highly reliable simulations on large lattices at 
low temperatures p3,E0]. These works indicate that, in 
fact, no QCP point exists below x p and that at perco- 
lation threshold the phase transition is characterized by 
classical exponents |^6|,[27| . 

Note that the theoretical studies of diluted spin sys- 
tems has attracted much attention some 30 years ago 
in the context of magnetism in diluted magnetic alloys 
p(J| p3[ . Most of these studies focused on the large S 
(classical) Heisenberg or Ising systems. Traditional view 



of the effect of local disorder on the spectrum of an or- 
dered 3D antiferromagnet is that at long wavelengths 
the form of the spectrum is not modified. The only ef- 
fects are the reduction of hydrodynamic parameters (spin 
stiffness, spin- wave velocity, etc. [B4J) and a weak damp- 
ing. Conventional arguments for this ability of the long- 
wavelength spectrum to withstand perturbations on a lo- 
cal scale often appeal to the Goldstone theorem |3J|, al- 
though its applicability to systems without translational 
invariance requires an additional assumption that the 
microscopic details are virtually averaged on short dis- 
tances. As a result the low-energy excitations are the 
weakly damped spin-waves which belong to the so called 
"infinite cluster" and they are well defined up to the per- 
colation threshold |35|] . This effective restoration of the 
translational invariance involves the spin-wave propaga- 
tion on randomly directed paths with some Euclidean 
distance L' which can be converted to a "true" distance 
L. Thus, the wave- vector preserves its meaning at long 
wavelengths pS. In 3D these arguments work very well 
and are assumed to demonstrate a "general principle" . 

There is growing evidence that in 2D, however, such 
a logic is not always valid. In the 2D case it was found 
by Harris and Kirkpatrick p0[ and more recently in Ref. 
p7|, using a perturbative (linear in x) approach, that 
the spin-wave self-energy at long wavelengths acquires a 
non-hydrodynamic contribution which explicitly violates 
the Lorentz invariance of the clean system (a feature an- 
ticipated by Chakravarty, et al. in Ref. Q). Recently, 
similar results have been obtained in RPA studies of the 
diluted 2D Hubbard model B. Some of these studies 
have concluded that D — 2 is the lower critical dimen- 
sion for this type of disorder [BO 37 1, implying an insta- 
bility of the long-range order to an infinitesimal doping 
in the Imry-Ma sense ]3q ]. However, as we mentioned 
above, MC results show that the order is preserved up 
to x — x p in contradiction with these conclusions. We 
will show that the conjecture of the instability is an arti- 
fact of a perturbative expansion and is avoided when the 
divergent series of diagrams is summed. However, the 
resulting modification of the excitation spectrum is very 
unusual and leads to a number of observable anomalies. 

Technically, our approach is similar to the one by 
Brenig and Kampf |L5| who have studied the problem of 
excitation spectrum in diluted 2D QAF using spin-wave 
and T-matrix formalism. However, while the authors of 
Ref. |1^] have noticed unusually broad peaks in the spec- 
trum, they assumed a "normal" 3D-type of the spectrum 
renormalization, that is, the softening of the sound ve- 
locity and a recovery of the spectrum at long wavelength. 
The derivative of the spin-wave velocity with x, obtained 
numerically in Ref. |15| using this assumption, is rather 
large d(c(x) / Co) / dx w —3 which has supported earlier 
experimental expectations of a QCP at x < x p . A recent 
work by two of us using the non-linear er-model allied 
to classical percolation theory |3$j] gave similar result. 
Another study in Ref. E(J used a generalizations of the 
cr-model with parameters modified according to MC data 



and suggested a simple renormalization of spin stiffness 
p s (x) and spin- wave velocity c(x) as the only effect of 
impurities. We will show that these results are not cor- 
rect because of the existence of localized spin excitations 
which are not taken into account in these works. 

In this work we study the problem of impurities in 2D 
QHAF within the linear spin-wave theory using the T- 
matrix approach combined with configurational average 
over the random positions of impurities. We solve the 
single impurity problem exactly. The spin-wave Green's 
function is evaluated by summing all multiple-scattering 
diagrams that involve single impurity. This approxima- 
tion gives results that go beyond simple linear expansion 
in x, although not all higher order contributions in x 
are taken into account. This approach is valid as long 
as single-impurity scattering is the dominant one. We 



recover the results of Rcfs. [B0U37H at k ^> 



that is, 



the spin- wave spectrum acquires a non-linear logarithmic 
contribution £k (<■<->) ex x k In \u>\ with an abnormal damp- 
ing Tk ex x k. This means that an effective spin- wave 
velocity c(x) is not well defined. However, we show that 
there is no instability of the system towards a disordered 
phase, as conjectured previously. The static properties 
such as staggered magnetization and Neel temperature 
do not possess anomalies in contrast with the dynamic 
properties. It is interesting to note that the spin- wave 
stiffness p s {x) — c(x) 2 x±(x) is also well defined since the 
anomalous terms in the transverse susceptibility x±( x ) 
and in c(x) cancel each other. 

We show that the diluted 2D AF provides an example 
of the system where the arguments for the spectrum to 
be "protected" at long wavelengths fail j|l],[42| . We have 
found that the spectrum of a 2D AF at long wavelengths 
is overdamped at arbitrary concentration of spinless im- 
purities. More explicitly, the spectrum ceases to contain 
a quasiparticle peak of any kind beyond a certain length- 
scale. The actual spin excitations instead of being de- 
scribed as ballistic may be interpreted as diffusive spin 
modes. The reason for that is the influence of scattering 
centers on the long-wavelength excitation which is not 
vanishing in 2D because of the small phase space. This 
leads to the absence of the effective self-averaging of the 
system to a translationally invariant medium with the 
renormalized parameters as it would be in 3D. Instead, 
the scattering leads to a new length scale £/a ~ e*/ ix be- 
yond which the influence of impurities on the spectrum 
is dominant. We associate this length scale with the lo- 
calization length of spin excitations. 

We show that the dynamical structure factor 5(k,w) 
for a -1 ;§> k ;§> £ _1 consists of three parts (we use units 
such that h = ks — !)'■ (i) a broadened quasiparticle 
peak at the energy u — Cok (1 + 2xln(fca)/7r), where 
Co = 2y2SJa is the bare spin- wave velocity, J is the 
antiferromagnetic exchange constant, with a width given 
by Tj, ~ xc^k] (ii) a non-Lorentian localization peak at 
ui — ujq ~ co^ _1 , (Hi) a flat background of states be- 
tween u> = cok and u> = ujq. Thus, besides the lack of 



the Lorentz invariance, for every k-state some weight is 
spread from the h igh energies ui ~ cofc to the low energies 
down to lj ~ ll>o |43[ . For fc < i? -1 the quasiparticle and 
localization peaks in <S(k, to) merge into a broad incoher- 
ent peak that disperses in momentum space. 

The anomalies in the dynamical structure factor are 
reflected in the magnon density of states, N(ui). In a 
clean 2D AF N(uj) oc w. With the doping N(u>) ac- 
quires a constant contribution from the localized states 
N(lu) oc lo + const ■ i at w > ujq, has a peak at u> rs loq 
of the height ~ 1/x, and vanishes as N oc l/(xln \w\) 2 
for cj ^C cjo- This behavior of N(u>) is reminiscent of the 
problem of localization of Dirac fermions in 2D <i-wave 
superconductors Q in the case of "strong" disorder (uni- 
tary scatterers). Another interesting similarity between 
that problem and impure 2D AF is that disorder may lead 
to very different physical consequences depending on its 
"class". As it was noted in Ref. |37J] and also in Ref. 
Pq| in another context, one obtains drastically different 
results if the spin of impurity is equal to the spin of the 
host material and only bond strengths around impurity 
(J) are modified. The renormalization of the spectrum in 
that case does not contain any anomalous terms, namely: 
Sk(w) ~ a; cofc and Tk ~ xcok 3 a 2 . According to the ter- 
minology of 2D Dirac fermions this problem falls into the 
class of a "weak" disorder. In the case of spinless impuri- 
ties the similarity to "strong" unitary scattering centers 
is evident since no spin degrees of freedom exist at the 
impurity site. 

From the density of states N(w) we calculate the mag- 
netic specific heat which for a clean 2D system at low 
temperatures is Cm(0, T) oc T 2 . We predict a strong de- 
viation from this behavior due to localized states. We 
find that specific heat acquires a quasi-linear correction 
SC M (x, T) = 0(x)T/( In 2 |T/w |+vr 2 /4) which is roughly 
oc xT at T ^ wq, j3(x) ~ 1/x. Observation of such a be- 
havior can provide a simple test of our theory. We remark 
that in our approach the contribution of the finite (decou- 
pled) clusters is not taken into account since the whole 
system is considered as a single, ordered, infinite cluster. 
However, finite clusters of the size L have a gap in their 
spectrum of order J/L and thus become important in 
the low-T region only at x close to percolation threshold 
where L can be large. Another source of similar high- 
energy corrections is from the resonant states (uj res ~ J) 
around impurities whose energy may go down with dop- 
ing |Q. At lower temperatures, T < y/J±J, where J± 
is the inter-plane exchange constant, the crossover to a 
3D behavior should be seen. Thus, for x not too close to 
x p we expect a large temperature window where the pre- 
dicted anomalous 2D behavior of the Cm in the infinite 
cluster is dominant and can be observed. 

We also consider the effects of small inter-plane cou- 
pling t 3 £> = J±/2J and small anisotropy gaps on our 
conclusions for dynamical properties of a strictly 2D 
isotropic AF we discussed above. It is evident that 
as long as these additional energy scales are small 
in comparison with J there will be an energy range 



1 ;§> uj/J » y/r (r = T e ff accumulating the total 
effect of the gaps and 3D coupling) in which the non- 
linearity of the spectrum and an abnormal damping of 
the 2D spin waves should be observable. A more deli- 
cate question is if the localization part of the spectrum 
and truly overdamped long-wavelength excitations can 
be seen in the presence of gaps or 3D coupling. The 
point is that the disorder induced scale loq ~ Je^^/ ix 
can be hindered by these additional terms which cuts off 
the log-singularity. Therefore a range of concentrations 
< x < x* ~ In (1/t) can be found where the long- 
wavelength quasiparticles are still well defined deep in 
the 3D region of the k-space (ka <C y/r) similar to the 
quasi-lD problem Q. For the LCO materials r ~ ICP 4 
gives x* ~ 0.1 — 0.2. Above the concentration x* (and at 
x < x p ) localization and overdamped peaks should be ob- 
servable since ujq > yfr and all the low-energy excitations 
become incoherent. Our order of magnitude estimation 
for the largest value of r which can allow such observa- 
tion (from the condition x* < x p ) is r ~ 0.01. Therefore, 
a rather high impurity concentration and small enough 
anisotropics and inter-planar coupling may be required to 
observe directly some of the dynamical effects we predict 
in this work. 

We calculate the static magnetic properties and find 
a quantitative agreement with both MC simulations and 
experimental data. We show that at T = the staggered 
magnetization (averaged over the magnetic sites |4q| ) is 
given by M(x,0) w S — A — Bx for x <C 1, the factor 
A = J2k w k ~ 0-2 stands for the contribution of the zero- 
point fluctuations of the spins, B ~ 0.21 is S-independent 
in our approach. We find that Tn(x)/Tn(0) ~ 1 — A s x 
for x <C 1 where A s = tt — 2/tt+B/(S—A) is a weak func- 
tion of S. This linear expansion result gives A± /% ~ 3.2 
and A 5 / 2 — 2.6 which work quite well up to a high value 
of x ~ 0.25. It is interesting that the linear expansion 
results point to x c (l/2) ~ 0.31 and x c (5/2) ~ 0.38, both 
below x p , which means that Tjv(x) versus x curve should 
be concave, in contrast with the 2D Ising magnets for 
which Tjv(x) is a more traditional convex curve p7j . Such 
an anomalous curvature of the ordering temperature has 
been also observed in many different magnetic systems 
composed of /-electron moments such as U and Ce Q . 
We show that in our approach for larger values of x 
Tjv(x) indeed bends inward and tends to saturate close 
to x p . We interpret this behavior as due to localization 
effects which tend to reduce the role of quantum fluctu- 
ations in the destruction of the long-range order. 

We have calculated the 2D magnetic correlation length, 
£(x, T), to describe the paramagnetic phase of the system 
above the Neel temperature. We used a modified spin- 
wave theory formalism by Takahashi |p0| and calculated 
£(x,T) numerically. Correlation length is suppressed in 
comparison with the pure case and also shows some de- 
viation from the simple e 27TPs ^ x ^ T behavior at larger x. 

This paper is organized as follows: we describe the 
model and introduce the formalism in Section O; in Sec- 
tion III, we present the results for the dynamic proper- 



ties; in Section IV the static properties and long range 
order is discussed; Section |yj contains our conclusions. A 
few Appendices are included with details of the calcu- 
lations. Some of the results presented here were briefly 
reported in our previous paper pjj. 



II. FORMALISM 



The systems discussed in this paper are modeled by 
the site-diluted quantum Heisenberg antiferromagnet: 



H = ^ J v Pi Pj Si • Sj , 



(i) 



where pi — 1 (0) if R^ site is occupied (unoccupied) by 
the spin S. We focus on the problem of tetragonal or 
square lattices with in-plane, J, and out-of-plane, Jj_, 
nearest-neighbor exchange constants, (ij) denotes sum- 
mation over bonds. In the systems of interest J ;» J± 
(for instance, in LCO J « 1500 K and J ± ~ 10~ 4 J). 



1. Spin-wave approximation. 

We begin with the Hamiltonian (pi) which is split into 
the pure host and impurity part 

H = H{)+ Himp = 2_^ Jij ^8 ' Sj — 2_^ 'h,S S; • S/ +< 5 , (2) 



m 



1,6 



where I runs over the impurity sites and S is a nearest- 
neighbor unity vector. Then, in the linear spin-wave ap- 
proximation, 



Sf 



S — a] a;, S) 



Sf = -S + btbj, S. 



2Sdi, S t 



>2Sb), SJ 



2Sa 



256 



3 ' 



(3) 



for the spins in A (i) and B (j) sublattices quadratic 
part of the pure host Hamiltonian H.q for the tetragonal 
lattice is given by: 



Ho = 4:SjJ2 



7o(a k ak + fo k fok) 

+ 7k(a k 6 -k + & -k a k) 



(4) 



where we use that in-plane and out-of-plane coordination 
numbers are z — 4 and z±_ = 2, respectively, and define 



7k = 7k + T7 k 



(5) 



with r = Jj_/2J, 7 k = (cosfc x + cosk y ) /2, and 7 k = 
cos k z . From now on the in-plane and out-of-plane mo- 
menta are in units of the correspondent inverse lattice 
constants. Impurity part of the Hamiltonian (0) on the 
tetragonal lattice is: 



Hfmp = S Yl Jl > S ^l ai + b \+S bl + & + a l h \+S + a ^+S 
1<£A,S L 



(6) 



with Jis = J (J_l) for 5 = e x ,e y (e z ). After Fourier 
transformation it is more convenient to write impurity 
Hamiltonian in the 2x2 matrix notations: 



U lmv = -45 J J^ e J ( k - k ') R 'i k ^ jk ,i k , , 



i,k,k' 



where 



A v = 



«k 



. 4 = [ 4, b 



with scattering potentials for I in the sublattice A: 



V A - 



7o 7k' 



7k 7k-k' 

and for I in the sublattice B: 

7k-k' 7k 



v B - 



7k' 7o 



(7) 



(8) 



(9) 



(10) 



The pure host Hamiltonian (0) is diagonalized using Bo- 
golyubov transformation: 



a k = M k a k + v k f3i k 
bl = Uk/? k + «ka_k 



with 



"k _ v k = 1 > 2u k«k = -7k/wk , 



(11) 



(12) 



/7o+^k . /7o-^k 

Uk = \/— ^ , «k = -sgn7 k1 ' 



2w k 



where bare spin-wave frequency is 



2wk 



"-"<< = \/7o - 7k 



(13) 



The problem can be reduced to the problem in 2D square 
lattice by letting r — ► in Eqs. (fC-(jl3|). In what follows 
all energies are expressed in the units of Oq = AS J. 

After the Bogolyubov transformation the Hamiltonian 
Eqs. fl),(H) is given by (in the units of Qq): 



Wo = J^cj k fa k a k + ^ k /3 k j , 



n mp = -J2 e I(k - k ' )R U k v k:k ,A' 

;,k,k' 

where two-component vectors are: 



A 



Ok 



A 



at, (3-i 



(14) 
(15) 

(16) 



and 2x2 scattering potential matrices V k k' are obtained 
from Eqs. (|),© using Eq. (fl|). For the sake of the 
further use of the T-matrix formalism it is convenient 
to decompose scattering potentials into the orthogonal 
components according to the symmetry with respect to 
the scattering site. The symmetry of the tetragonal lat- 
tice is D^hi which is a group of order 16 and has 10 ir- 
reducible representations. Since the impurity potentials 
Eqs. (^,0) connect only nearest-neighbor sites, only five 
components of the scattering potentials in the irreducible 
representations of D^h are nonzero. They correspond to 
the irreducible representations A\ g , B\ g , B? g , and E u . 
These nonzero components are the s-wave, in-plane p x -, 
p y -, and (i-waves, and out-of-plane p z -w&ve (for details 
see Appendix [A|). 

Thus, the scattering potential for the impurity in the 
sublattice A: 



V A - 
nc.k' — 



EK$ 



(17) 



where scattering channels are H = s,p x ,Py,d,p z . In each 
channel the scattering potentials can be written as a di- 
rect product of the column and row vectors. The s-wave 
part 



problem further by neglecting the s 1 - and p z components 
in the above equations. Moreover, the solution for the 
2D problem can be applied directly to the quasi-2D case 
since the formal expressions are identical in both cases. 
The only important difference concerns the logarithmi- 
cally divergent terms which in quasi-2D system acquire 
a low-energy cut-off provided by the implicit dependence 
of the scattering potentials (|19|)-(|22|) on r through Wk- 
This simply means that for the quasi-2D case in the 
limit t«1 one can restrict oneself by considering purely 
2D scattering including the 3-dimensionality only on the 
level of the spin-wave dispersion in certain terms. Thus, 
in the following we use 



V, 



k,k' 



Sk) ® (Sk' 



with (s k | = cj k ["k, -«k] 



(22) 



The rest of this section is devoted to the 2D limit of the 
problem and, unless specified otherwise, we use 



7k = 7k = (cos k x + cos k y ) /2 
7o = 1 , w k = V 1 - 7k • 



(23) 



VQy = W)®{sv\+T\si)®{si,\ , 

where (s k | = [W + « k 7k, w k + Uk7k] , (18) 



and 



[u k + «k7k > V k + Uk7k] : 



the in-plane p-wave part 

Vk^ w = bk (y) )®(Pk^i. 

where (p^. (y) \ = smk x{y - ) [v k , Uk]/V2, (19) 
the d-wave part 

VjJ£ = |dk>®<dk'|, 

where (dk| = 7k [ u k, u k ], (20) 

with 7 k = (cosk x — cosk y )/2, 

and the out-of-plane p z -wave contribution 



V, 



where (p£| = sinfc z [v^, u k ]/V2. 



(21) 



B,fi 



'}A,H I 



For the impurity in B sublattice V k '£, = V k k /(i< «-> v). 

In what follows we consider the 2D (r = 0) or quasi- 
2D (r€l) limit of the problem. It can be shown that 
the contribution of the out-of-plane terms in s-wave and 
p z -wave scattering potentials which explicitly depend on 
r, as well as the one of the majority of the r-dependent 
terms originating from the quasi-2D form of u^, «k> and 
Wk (p"2|), ( p^| ) is negligible in the quasi-2D case (~ O(r), 
see Appendix 0) . It allows one to simplify the scattering 



2. T-matrix. Single-impurity scattering. 

We are interested in the Green's function of the Hamil- 
tonian (114) modified by random impurity potentials (fi~q). 
The Green's function is a 2 x 2 matrix defined in a stan- 
dard way: 

Gi\t) = -i<T[ok(t)4(0)]), 
Gi\t) = -»(r[a k (t)0_ k (O)]), 
Gi 1 (t)=-i(T{pl k (t)ai(0)]), (24) 

G k 2 («) = -z(T[/3 k (f)/3 k (0)]>, 

where brackets also imply a configurational average over 
the impurity sites. 

The T-matrix equation for the Hamiltonian (Hq) is 
given by 



^CM 



-V, 



k',k' 



T.^G»T^) 



(25) 



with ^ = A(B), partial waves are restricted to in-plane 
\x = s, p a , d harmonics according to the above discussion, 
and Gq(w) is the 2x2 bare Green's function: 



GiH-u) 



-(0,11/ 
q \~/ ~q 

q ( w ) = G°- a V)=<>. 



w).ia 



1 



W — Un + i0 



(26) 



The diagrammatic equivalent of the Eq. (B5h is shown in 
Fig . nVa). The T-matrix equations (Eq) with potentials 
(0)-(B2h can be readily solved: 



~-iA,fi 






- V B ^7 



M, 



(27) 



where the frequency dependent parts are given by 

1 (l+q;)pH 

S[ } lu + 1- W (1 + w)p(w) ' 

r , s 2 

pW l + w + (l-w)(w*p(w)- Pd (w))' 

T d (w) = - 

with 



(28) 



1 + (1 - oj)p d (uj) ' 



Pel 



H = E^r- <-'»> 



W 2 — U 2 



We note here that the second term in s-wave T s (uj) is 
proportional to p(oj) at u> <C 1, where the latter appears 
naturally from the summation in Fig. 0(a) as a result 
of combination of (7 q ' u (u;) and G q ,22 (w) in the internal 
part of the diagrams. When the summation over p in 
Eq. ( P9[) is restricted to 2D p(u) is a logarithmic func- 
tion at low energies. In the following we show that this 
contribution to the s-wave scattering is solely responsi- 
ble for all the anomalies in the spectrum of a 2D AF. 
Interestingly, similar logarithmic term in the self-energy 
of the 2D Dirac fermions in the problem of disorder in 
d-wave superconductors requires a summation of the spe- 
cific subset of diagrams [p2|. In our case, while one needs 
to sum infinite series of diagrams, no special selection or 
inclusion of the multiple-impurity scattering processes is 
necessary. Since the single-particle density of states and 
sensitivity of the results to the type of disorder in both 
problems are similar, establishing of the detailed corre- 
spondence between these two problems is an important 
question. Integrals in Eq. (E3) can be taken analyti- 
cally and, in the case of 2D, are expressed through the 
complete elliptic integrals [15[ (see Appendix O . 

The first term in the s-wave scattering Eq. ( |28| ) rep- 
resents a singular zero-frequency mode which is indepen- 
dent of the dimension of the problem and originates from 
the oscillations of the fictitious spin degrees of freedom 
at the impurity site which are decoupled from the AF 
matrix. Roughly speaking, when the spins are quantized 
as in Eq. (0) and S' at the impurity site is set to zero 
there is still a a,Q left from Sq. Thus, in the spin- wave 
approximation, it gives rise to the s-wave zero-frequency 
mode. This problem has been noticed since the earliest 
works on the diluted magnets which have used the spin- 
wave theor yjplj and also more recently in the context of 
diluted AF fl7j, |15| , [^ , [45|]53| l (for an extensive discussion see 
Ref. p7[ ). Since these states are unphysical and are unre- 
lated to the low-energy physics of the AF they have to be 
projected out. One of the projection schemes involves a 
non-Hermitian potential which was designed to preserve 
the simplest factorized form of the s-wave scattering po- 
tential pTfl. We use another, physically more transparent 



scheme which introduces a fictitious mag netic fields at 
the impurity sites (similar to Refs. |15| . [l5| ): 



AH = H Z J2 a]ai 



(30) 



i,k,k' 



,i(k-k')Bj it av'' s 



41 Ay L '» A, , 



where corrections to the s-wave scattering potential are: 



Av£& = |As k ) <g> (Aa k , 



with (As k | = [u k , Vk] , 



(31) 



^Hc ic' = ^kit'l" *""* v }- Evidently, p and d waves 
are not affected by the projection. Within our approach 
after some algebra in the limit H z — ► 00 one obtains a 
modified T-matrix solution (for the case of an arbitrary 
H z see Appendix y ) : 



t4>) = to s M + Af£>) 



T, 



k.k 

B 

k,k 



(32) 



^)=V k 3 £T s (-u) + Affi( u l 



where V kk / is given, as before, by Eq. (|l^) and the 
frequency dependent part is now free from the zero- 
frequency pole 



r a (u) 



{l+u))p(u) 

1 — Lj(l + (jj)p(ljj) 



(33) 



Comparing this expression with Eq. ( p8f ) one may note 
that the "physical" term is left unchanged after the pro- 
jection. Additional terms in the solution ( p2[ ) are also 
regular: 

A ?£k'M = -w|As k ) ® (As k ,| (34) 

+ |sk)®(Ask<| + |Ask)®(sk'| , 



with |s k ) from Eq. 



and |As k ) from Eq. Q, 



i k jj(u) = T w £,(—u>){u <-» v} as before. Thus, the pro- 
jection ( J30| ) allows one to remove the unphysical diver- 
gency at ui — which would otherwise affect the true 
low-energy physics of the problem. 



3. Green's function. 

The averaging over random distribution of impuri- 
ties readily transforms T-matrix into the spin-wave self- 
energies: 



Sk(w) = y^S M:k (cj) 



with /z-wave contributions 



(35) 



(36) 



For the 2D case the contribution of the partial waves to 
the self-energies are: 



S a ,k(w) = XLU k 



i 7 |,\ r,H + r,H 
7k i J 2 

^k o ~\ r a (w) - r 8 (-w) 

-Wk 



(37) 



S p ,k(w) = iwk 



2 
2 



- \ 2 

2k 



2 

1 
-1 



(38) 



S<j,k(w) = iwk 



2k 



1 - 7k \ r p (o;)+rp(-a;) 

-7k i y 2 

-Wk \ r p (o;) - T P (-uj) 
uj k J 2 

1 - 7k ^ T d (Lj)+T d {-u}) 

-7k 1 



-w k \ T d {uj) -r d {-tu) 
cj k 



(39) 



It is interesting to observe that "on shell" (at u> = w k ) 
"projected" T k , k ,(w) from Eqs. @-@) and "non- 
projected" expressions Eq. (p7|), (|2§|) yield identical 
£g,k(wk)- 
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FIG. 1. (a) T-matrix single-impurity scattering series, (b) 
Dyson-Belyaev diagram series for the diagonal, G , and 
off-diagonal, G 12 , Green's functions. Self-energies E 11 (cir- 
cle) and E 12 (square) are the configurational averages of T 11 
and T 12 components of the T-matrix, respectively. 



Summation of the Dyson-Belyaev diagrammatic series 
for the Green's functions shown in Fig. pl(b) with self- 
energies defined in Eqs. (p5|)-(^9|) gives 



G k (w) 



-u> — a»k- S k 2 (w) 



u,-SW) I !J! " 



(W - Wk 



£ k ») (-w - ^k - s k 2 H) - (££»)' 



where £ k 2 (w) = S k 1 (— uj). A detailed consideration of 
the properties of spectral functions 



4» 



i 



-ImG^M , 



(41) 



will be given in the next section. 

We investigate the neutron scattering dynamical struc- 
ture factor, 5(k, w): 



5°"(k,w)= / 

J — ! 



dte &rt <S£(t)5^ k (0)) , 



(42) 



which is directly related to the spin Green's functions. 
The standard derivation of the single-magnon contribu- 
tion to the transverse component of the dynamical struc- 
ture factor iS H (k, u>) at T — gives 

S + ~ (k, u) = n S (?i k + u k ) 2 



KH+ A k 2 H + 2A k 2 H] 



(43) 



where the kinematic (w-independent) form-factor (w k + 
v k ) 2 = (1 — 7 k )/w k is proportional to k close to the "nu- 
clear" reciprocal lattice point K = and is ~ 1/fc close 
to the "magnetic" Q = (tt, it) point. It thus enhances 
the signal close to the AF ordering vector and suppresses 
it close to the zone center. Note that the diagonal parts 
of the Green's function are symmetric and off-diagonal 
parts are asymmetric with respect to the transformation 
k -» k + Q (since G 12 ~ E 12 and S k 2 +Q = -£ k 2 ). There- 
fore, the sum of the spectral functions in the bracket in 
Eq. (J43|) is, generally speaking, different in the magnetic 
and nuclear parts of the Brillouin zone. At T > above 
expression (43) is modified by the factor [1 + tib(w)], 
where ns(w) = [e u / T — 1] _1 is the Bose distribution func- 
tion. 

The density of states associated with magnetic excita- 
tions is straightforwardly related to the magnon Green's 
function (f40|) and is given by 



tf( W ) = £[^"(" 



^ 2 M] 



(44) 



The magnetic specific heat is then given by 
1 



C M {T) 



rp2 



dLjN(uj)uj 2 [n B (uj) 2 +n B (u;)] , (45) 



where w and T are in the units of f2o = 4S 1 J. 

The static properties of the system, such as staggered 
magnetization in the ordered phase, Neel temperature, 
and 2D correlation length in the paramagnetic phase, are 
calculated from the spin- wave expression of the averaged 
on-site magnetic moment: 



1 



\(Sl)\=S-Wx 

2 k v^ 






(46) 



-7k (44) 



where bosonic averages can be expressed through the 
spectral functions (|4l| ) as: 



(ajc"k) 



du> tib(w) 4 Hk (w) 



(at(3l)= J cb;n B (u)A% k (u;), 



(47) 



which implicitly depend on impurity concentration x, in- 
dex R stands for retarded. 

In the ordered phase these expressions (|46|), ([47]) pro- 
vide us with the concentration and temperature depen- 
dence of the averaged staggered magnetization M(x,T). 
The same expressions with the condition (S z )(x,T) — 
define the mean-field equation on the Neel temperature 
as a function of a;. In both cases, when T ^ the 3D 
form of the spin-wave dispersion is to be used in Eq. 
@). In the paramagnetic phase (T > T N ) Eq. Q ) 
should be modified by 7k — > 777k and Wk — > yl — ?? 2 7k- 
Then, in the framework of the modified spin-wave the- 
ory I5CJ1 , equation (S*)(x,T,ri) = is a constraint which 
represents a self-consistent equation on the gap y 1 — if- . 
This, in turn, defines the 2D correlation length £2.0 as a 
function of 2: and T. 



III. DYNAMIC AND THERMODYNAMIC 
PROPERTIES 



In this Section we consider in detail the structure of the 

fectral functions of the Green's function Eq. (EO), Fig. 
b), with self-energies given by Eqs. (p5|)-(|3^). We cal- 
culate the dynamical structure factor <S(k,w), spin- wave 
density of states N(u>), and low-T magnetic specific heat 
Cm(T). We consider the long- wavelength, low-energy 
limit of the problem and obtain analytical results for the 
low-energy <S(k,w) and N(u>) and the low-temperature 
Cm(T). We recall here that all wave- vectors are in units 
of inverse lattice spacing 1/a and all energies are in units 
of Q = <±SJ. 

We consider the low-energy form of the Green's func- 
tions first. At low energies CL>,Wk <C 1 self-energies (p7|)- 
( p9| ) are given by: 

X^M = iwk^w) + 2 - tt/2] - xl) + C(w k w 2 p 3 ) , 
£ k 2 H = zw k 7 k \p{u) + tt/2] + 0(«kwV) , (48) 

with p(u) ~ (2/7r) In |w/4| - i , 

which includes contributions from s- and p-wave scat- 
tering, d-wave part is of higher order E^ ~ O(oj^), 
Wk — kj\[2. The importance of the projection of the 
unphysical states can be demonstrated one more time by 
comparison of the above expressions with the "unpro- 
jected" (H z = 0) form of the self-energy: 

E£ x (w) = xl)^[p(lj) - tt/2] + xluI/lu + 0{u^ 2 p z ) , (49) 



which possesses an ui — singularity. Noteworthy, the 
"physical" part of the expression containing the loga- 
rithm is not related to the unphysical states and stays 
intact under the projection. As it was noted before "on- 
shell", u) = cjk, the self-energies (E§) and p9) coincide 
|p7| . The off-diagonal S k 2 (w) is the same in both cases. It 
is also useful to note that the first-order Born approxima- 
tion to the scattering problem would give a very different 
result 



rill, Bom/ \ o 



s 12,Borr l(w) = Q ^ ^ 



with imaginary part of the self-energy being ~ 0{xkuj 2 ). 
One can see that along with the "normal" softening 
and weak damping the full T-matrix consideration gives 
non-linear dispersion term and damping |7k|/wk — x 
which is only parametrically small with respect to the 
bare spectrum. A perturbative "on-shell" pole equation 
gives 

^k + «7k = w k + Sk 1 K) , 

W k =Wkfl- x{ir/2 - 1) + ^ In |w k /4| J , (51) 

7k = -xuj k , 

which already shows that the spin-wave velocity in the 
effective medium is not well defined since the bracket in 
Eq. ( pl| ) depends on k. Moreover, the renormalization of 
the real part of the spectrum is dominated by the In |w| 
term at low frequencies and the bracket vanishes at some 
wave- vector 



exp(7r/2x) 



(52) 



Because of that one can naively suggest a vanishing of 
the spectrum p7| and an instability of the ground state 
towards some new phase. Such an instability is, of course, 
just a signature of the breakdown of the perturbation 
theory. One has to sum up all the "dangerous" terms 
using Belyaev- Dyson equation Fig. Jj(b) and Eq. ( |40| ) 
and analyze the spectral functions (|ll ) . 

The low-energy, long-wavelength form of the Green's 
functions Eq. (M) with self-energies from Eq. (|48|) is 



G k 1 M=G k 2 (-o;) 

G k 2 H^ 



D + w k ( 1 + x[p(u>) + 2 - tt/2] 

Cj 2 - ujI (l + x[2p(uj) + 4 - tt] 
S7kW k [p(w) + 7T/2] 



Co 2 - ul 1 + x[2p(uj) + 4 - tt] 



(53) 



where w = w(l + x) and p(u>) is defined in (J48). 

The diagonal spectral function in the same limit can 
be then written as: 



I xw k (w + a;k] 



4c 1 M 



k (u> 2 - u>la(u>)y + (2xu)lY 



(54) 



where we make use of imaginary part of the self-energies 
being ImE^(w) ~ -iwk, and introduce a "stretching 
factor" 



o(w) = 1 + x(- In |w/4| + 4 - 7r) 



(55) 



The energy at which this factor vanishes defines the 
disorder-induced energy scale 



(56) 



Wo ~exp(- — 



below which the spectrum is overdamped. 

A more detailed analysis of Eq. (|54|) gives the follow- 
ing picture. At the wave-vectors much larger than loq 
(wit ^$> Wo)> that is at the wavelengths shorter than a 
characteristic length £ ~ e _7r / 4:r , the spectral function 
has three distinct regions in lo. First, is a vicinity of a 
quasiparticle peak, lo ss lo^: 



AlH") 



2w k 



2xw£ 



- (^ - ulY + (2xuiy 



(57) 



where the spectrum has a regular Lorentzian form with 
the pole at u>k and width 7k given by the perturbative 
result Eq. ( |5l] ) ■ Second, the intermediate range of ener- 
gies, uiq < lo <C Wk, where the "stretching factor" is not 
too close to zero: 



4c 1 M 



1 X 



1 



1 X 



IT LOk a(w) 2 + Ax 2 IT LUk 



const 



(58) 



one can approximate a(ui) by a constant since its depen- 
dence on lo is weak in this range. One can see that the 
spectral function in this region is independent of u and 
corresponds to an almost flat, shallow (~ x) background 
of states. Third, the vicinity of a "localization peak" , 
to kj luq: 



A" 



M 



Ai x M 



i i 

47T IWk 
7T 1 



at oj 



1 



16 w k a; In 



= wo , (59) 

at lo <C cjo , 



where the spectral function rises sharply from the shal- 
low background states ~ x (pq ) to a peak of the height 
~ 1/x and then vanishes in a singular fashion as lo ap- 
proaches zero. Note that this peak is non-Lorentian and 
its position (lo = loo) is independent of the value of k. 

Thus, besides the lack of the Lorentz invariance of the 
quasiparticle part of the spectrum Eq. (|5l]), every k- 
mode redistributes some of its weight from the energy 
~ cjk to a flat background of states between lo^ and loq 
and to a peak at lo = loq. Such a behavior is similar to 
the other problems of linearly dispersive excitations in 
the presence of disorder in two dimensions and should be 
interpreted as the signature of localization Q,Q . Then 
the characteristic length 



exp 



UxJ 



(60) 



is to be understood as a localization length of the spin- 
waves in our problem. 

The truly intriguing question is what is happening 
at the wavelengths of the order of I and beyond. In 
our approach for k < £~ x the quasiparticle and local- 
ization peaks merge into a broad incoherent peak that 
disperses in the momentum space. One can see that at 
k ~ loq < £ factor a(co) is negative and the "pole" in 
Eq. ( p4| ) becomes pure imaginary. However, since a(co) is 
lo dependent this peak is non-Lorentzian and thus can not 
be associated with the "simple" diffusive mode. Thus, we 
observe an overdamped, non-Lorentzian diffuse-like exci- 
tation with a characteristic width of the order of u>k and 
peak position roughly at U) < Wk- We have to remark here 
that the nature of the states at the wavelength above the 
localization length might be beyond the ability of our ap- 
proach and the proper description of them may require 
a different, non-perturbative type of study. 

Thus, the structure of the spectral function we dis- 
cuss above demonstrates an unusual, non-hydrodynamic 
type of behavior of the spin-excitation spectrum of a 
diluted 2D AF. The strong influence of disorder in the 
low-energy excitations in 2D results in the failure of the 
averaging procedure, which effectively restores transla- 
tional invariance, to recover the long-wavelength exci- 
tation spectrum of this effective medium. Already at 
the energies much larger than the disorder-induced scale 
lo ~ k ^$> loq one finds a departure from the hydrodynam- 
ics: while the "quasiparticle" excitation can be found, it 
does not disperse linearly with k and its damping is nei- 
ther hydrodynamic nor quasiparticle-like. More impor- 
tantly, above the characteristic wavelength £ no hydro- 
dynamic description of excitations is possible. The low- 
frequency modes do exist in some form but they cannot 
be classified in terms of an effective wave- vector and thus 
the long-wavelength propagation is entirely diffusive. 

In addition, the spectra at k 3> loq are not exhausted 
by the quasiparticle peak. They also consist of the back- 
ground of localized states and a localization peak de- 
scribed in Eqs. ©, @. 

The spectral function A\^(lo) obtained from the "full" 
expressions for the Green's function (Boh and self-energies 
(p5|)-(|39|) without taking the low-energy limit is shown 
in Figs. Eh3 for a number of wave-vectors along the 
(1,1) direction of the Brillouin zone for a representa- 
tive value of impurity concentration x = 0.1. The pur- 
pose of these pictures is to demonstrate the features we 
have discussed using the long- wavelength form of A\^ (lo) . 
The amplitude of each A^(lo) curve is normalized to 
fit the picture and therefore the relative heights of the 
curves bear no meaning. These figures also show the 
bare spin-wave energy (dashed-dotted line) with arrows 
pointing down showing the positions of "unperturbed" 
delta-function peaks. Dashed line corresponds to the per- 
turbative renormalized spin- wave dispersion, Eq. (|5l|), 
while arrows pointing up show the actual positions of 
the peaks for selected wave- vectors. The figures show 
the spectral function within the different ranges of k rel- 



ative to wo, k ^> ujq, k > loq, and k < l uq, respectively. 
The latter can be calculated using Eq. ( |56| ) which gives 
cj (a; = 0.1)-10- 3 . 



o.o 



o.oo 



0.2 



0.4 -i 



0.6 



0.8 



1.0 J 















x=0.1 


IS 












A^(co) - 


- 


N;x 


^x\ 








k » co 


] 




N 


s 


^T — 




- 








r- 








1 








\ v. 

N 

t \ 


\J> 




1 




i 






\ 


i >4 V 

i 


f 





0.0 



0.2 



0.4 . _ T 0.6 



o.,s 



1.0 



FIG. 2. The spectral function A k (uj) for the wave- vectors 
k = 0.1,0.3,0.5,0.7, and 1.0, all > uj along the (1,1) di- 
rection, k is in units 1/a. Dashed-dotted line is the bare 
spin-wave energy, arrows pointing down are the positions of 
original delta-function peaks. Dashed line is the renormal- 
ized spin- wave dispersion Eq. (plf), arrows pointing up show 
the actual positions of the peaks for selected wave-vectors. 
Aj. 1 (lu) for each k is normalized to fit the picture. 



Fig. g shows the spectral function A^ 1 (u>) for the wave- 
vectors k = 0.1,0.3,0.5,0.7, and 1.0 along the (1,1) di- 
rection (k is in units of 1/a, so that the corner of the Bril- 
louin zone is (tt,tt)). One can see that the quasiparticle 
peak follows the renormalized spin-wave dispersion ( |5l| ) 
at low k very closely. At higher values of k the higher en- 
ergy sub-band develops and the spectrum evolves into a 
"camel" -like structure discussed extensively in Ref. [ [l5| . 
The origin of this high-energy structure is in the presence 
of the high-energy resonance state (w res ~ J) around im- 
purity pa which is unrelated to the low-energy physics of 
the system. Since our low-energy consideration does not 
take this high-energy feature into account the position of 
the lower peak in this structure deviates from the long- 
wavelength dispersion d5l]) at larger k. The low-energy 
localization peak and background are already noticeable 
in Fig. g despite the high-energy scale. 

Fig. H shows the spectral function A^ (m) for the wave- 
vectors k = 0.005, 0.01, and 0.02 with the smallest wave- 
vector being of the order of ojq. One can clearly see the 
features we have discussed in Eqs. (57)-(p9|): the broad- 
ened quasiparticle peak, the localization peak, and the 
states between them. The quasiparticle peak continue 
to follow the renormalized spin- wave dispersion (fill) . As 
the k decreases all the mentioned structures merge. 
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FIG. 3. The spectral function A^ (o>) for the wave- vectors 
k = 0.005, 0.01, and 0.02 along the (1, 1) direction, k = 0.005 
is of the order of wo- Dashed-dotted line, dashed line, and 
arrows are as in Fig. H. Al^(co) for each k is normalized to fit 
the picture. 
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FIG. 4. The spectral function A^ (u>) for the wave- vectors 
k = 0.0001, 0.0005, and 0.001, all < u , along the (1, 1) direc- 
tion. Dashed-dotted line, dashed line, and arrows are as in 
Fig. ti. A^(cu) for each k is normalized to fit the picture. 

Fig. H shows the spectral function A^ 1 (u>) for the wave- 
vectors k = 0.0001,0.0005, and 0.001, all are smaller 
than loq. As we discussed above, the quasiparticle 
and localization peaks merge and give a broad, over- 
damped, non-Lorentzian diffusive peak. In other words, 
one may not represent the Green's function in this re- 
gion as a sum of coherent and incoherent contributions 
G^ h (uj) + G^ coIi (lj), it seems that only the second part 
survives. The peak position deviates from the pertur- 
bative renormalized spin- wave dispersion (51) and thus 
indicate the region where the perturbation theory breaks 
down. 

The off-diagonal spectral function A^(lu) should pos- 



10 



sess features similar to the one of the diagonal spec- 
tral function. The low-energy, long-wavelength form of 
A^(lo) is given by 



A 12 



k M * - 

7T 



1 xikUk(ujlb(x) - Lb 2 ) 



,7,2 



\a(u)Y 



(2xloI) 



(61) 



where b(x) = [l - 2x(n - 2)]. Note that A£{u) is not 
a positively defined function, it changes sign as a func- 
tion of lo at u> = LO^^J b(x) / {1 + x). Another important 
difference from A^(lo) is that A]^(lo) is odd under the 
transformation k > k + Q and thus has opposite sign in 
the first and second magnetic Brillouin zones. 

A detailed analysis of A^(lo) in different regions of 
<x>k and lo shows that in the vicinity of a quasiparticle 
peak A^(lo) has an additional smallness of order x in 
comparison with A^ (lo) , but it is of the same order in the 
"intermediate" (lo < w^) and low-energy regions where 
it can be approximated as 



At 2 1 



co) 



1 X7k 



n LOk a(io) 2 + Ax 2 



(62) 



with the behavior above, at, and below the localization 
peak identical to the one of A^-(uj), Eqs. (|5q), (p9|). 
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FIG. 5. The spectral function A^(uj) for the wave- vectors 
k = 0.1,0.3,0.5,0.7, and 1.0, all > lo along the (1, 1) direc- 
tion. AI?(lo) for each k is normalized to fit the picture. 

Fig. H gives an example of the structure of the off- 
diagonal spectral function A^(lo) obtained from Eqs. 
( [40| ) and (13a)- (39) without taking the low-energy limit 
for the wave-vectors k = 0.1,0.3,0.5,0.7, and 1.0. The 
features discussed in the preceding paragraphs such as 
change of the sign, low-energy localization peak, and 
background states are clearly seen in this spectral func- 
tion. 

The transverse component of the neutron-scattering 
dynamical structure factor S^ (k, lo) is directly related 
to the linear combination of the magnon spectral func- 
tions A£(w), A 2 ^(u)(=A^(-u>)), and A^(lo) as given 



by Eq. ([43]). It therefore must contain all the features of 
the spectral functions we discuss here. Fig. |6j shows an 
example of our result for iS" 1 (k, lo) v.s. lo at k = 0.1, 
that is in the "nuclear" Brillouin zone, for x = 0.05. 
Long dashed arrow shows the initial position of delta- 
functional peak. Since ±wo are very small in this case 
the localization peak is seen as a single spike at lo = 0, 
but the flat background of states is clearly visible below 
the quasiparticle peak. 



x=0.05 

k = 0.1 (1,1)" 




0.00 



0.05 /cr 0.10 
(it/zSJ 



0.15 



0.20 



FIG. 6. Transverse com- 

ponent of the neutron-scattering dynamical structure factor 
<S H (k,oj), k = 0.1, x — 0.05. Long dashed arrow shows the 
initial position of delta-functional peak. 

However, the actual observation of anomalous features 
of the spectra can be complicated because of two reasons. 
First, the structure factor contains a kinematic form- 
factor which enhances the spectral function combination 
by ~ 2/wk close to k = Q and suppresses it by ~ oo^/2 
close to k = 0. Second, as we show below, the sum of the 
spectral functions entering <S H (k, lo) is "less anomalous" 
close to k — Q than at k — > 0. Namely, the quasipar- 
ticle part of the spectrum is abnormally broadened and 
disperses nonlinearly for both k — > and k — > Q, while 
the low-energy localization features are suppressed in the 
vicinity of Q due to cancellation between the diagonal 
and off-diagonal contributions. 

One can show explicitly using the low-energy, long- 
wavelength limit of the sum of spectral functions 



[A£(w) + A1 2 (lo) + 2Al 2 {u)] given by 



A^(co) 



* E 

a/3=i,: 

Co 2 (l- 



.a/3 



O) 



2xL0k 



(63) 



7k) + w£(l + 7k) - 2x7kw£(7r - 2) 



i(lo)) + (2xto 2 ) 

that aside from the kinematic form-factor the dynamical 
structure factor should be different in the first (k — ► 0) 
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Al 



M 



Axlo^ 



t (&> - ula{u)y + (2xu£Y 



and the second (k — ► Q) magnetic Brillouin zones 



^H«z 



Axuj^uj 2 



* {Co 2 - w£o(«)) 2 + (2x W 2) 2 



(64) 



(65) 



due to the asymmetry of A£. 2 (w) to k — > k + Q. One 

can see that around the quasiparticle peak lo ~ lo^ these 
expressions are identical and are simply equal to the di- 
agonal spectral function Eq. (pTj), but at lower energies 
for k ~ Q the localization features are suppressed by the 
factor of lo 2 . 

This asymmetry is demonstrated in Fig. M which shows 
the intensity map of S^ (k, lo) -a;k/(l — 7k) = TtSA^ioS), 
that is the structure factor divided by the kinematic 
form- factor, in the k— lo plane across the Brillouin zone in 
the (1, 1) direction from k = to k = (tt,tt), for x = 0.25. 
The higher intensity corresponds to the higher value of 
the function. One can clearly see all the features of the 
spectrum described in this Section: the resonance and its 
splitting from the dispersive mode at high energies, the 
low-energy damped spin-wave mode in both the center 
and the corner of the Brillouin zone, and the asymmet- 
ric background of localized states with the low-energy 
peak at the bottom. The nonlinearity of the quasiparti- 
cle mode also seem to be quite visible though the actual 
detection of it or of the abnormal fc-dependence of the 
damping can be a challenging experimental problem. 



Evidently, the low-energy localized states should strongly 
affect N(ui) and one readily finds the anomalous correc- 
tions already on the level of perturbative analysis of the 
Green's function. If one uses the full T-matrix form of 
the self-energy but expands the Green's function in x: 



u k 



M 






•g£»e 



(u;)G?' u (uO 



one immediately gets a constant correction 



N(u) = -lo + xC + 0{xcu In \u>\) , 

7T 



(67) 



(68) 



which also implies a finite density of states at w = 0. 
A more sensible result can be obtained without using x- 
expansion from the long-wavelength expression for the 
spectral function A^-(co) (p4|): 



N(w) = -uj + xCi/[a(uj) 2 + Ax 2 ] + 0(a;u;ln|w|) , (69) 



where a(ui) is the same "stretching factor" Eq. ( |55| ) 
we used in Eqs. fj), ©, ©-©. At lo > lo q 
a(uj) pa const and we are back to the previous result given 
by x-expansion perturbation theory (|68|). At lu ps ujq den- 
sity of states has a peak of the height ~ l/x whose origin 
is evident: the low-energy non-dispersive localized states 
altogether contribute to it. At w < ujq the density of 
states vanishes as N oc l/(cc In |cjj|) 2 . Such a strong de- 
pendence of the result on the degree of approximation is 
reminiscent to the dispute over N(lo) for the certain types 
of disorder in 2D systems with linear excitation spectrum 
p4,p2| where different approaches result in drastically 
different answers for the low-energy part of the density 
of states. 



FIG. 7. Fhe intensity map of wSA k (uj), in the k — lo plane 



for k from (0,0) to (-7T, 7r) in the (1,1) 
u> = to u) — 1 for x = 0.25. 



direction and from 



The density of states of spin-excitations can be easily 
calculated using Eqs. (44) and (f40|). We recall that for 
the pure 2D system with linear spectrum of excitations 
low-energy density of states is a linear function of lo and 
in our case 



N(oj) = -lo . 
n 



(66) 




0.25 



FIG. 8. Density of states N(u>) v.s. lo for x — (pure 
system, dashed curve), x = 0.1, and x = 0.2 (dotted and 
solid curves). Dotted curves are the long- wavelength result 
Eq. (|(3|) with C\ — 4/tt' a/2 . Solid curves are the result of 
numerical integration using the full Green's function mu). 
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Our Fig. H shows the results for the density of states for 
x = (pure system, dashed curve), x = 0.1, and x = 0.2. 
The dotted curves show N(ui) given by Eq. (69) with 
C\ = A/tt^/ 2 which is obtained from a long-wavelength 
expression for the spectral function (|3) . The solid curves 
are the result of numerical integration using the "full" 
Green's function (40). While the overall agreement of 



wavelength expression for the density of states we obtain 
for such a correction: 



these curves is very good there is a significant discrepancy 
at low energies which has the following origin. In the 
long-wavelength limit we regarded the localization peak 
at lo = ujq as non-dispersive, whereas at larger k, close 
to the magnetic Brillouin zone boundary, it does disperse 
down to u ~ uj 2 ~ e~*l 2x <C u>q. This can be noticed 
in our Fig. p] as well. As a result of such a dispersion 
the peak in the density of states at ujq is spread to lower 
energies. Technically, there is a term in denominator 



of the Green's function ~ x p(u>) 



negligible at low 



k, which leads to such a behavior. Since this term is 
of the order of x 2 and our approach does not take into 
account all such terms we have no certainty on whether 
it is a spurious feature or not. As we show below this 
discrepancy does not affects any of our conclusions. 

In this context it is interesting to note that the con- 
stant term in the density of states, which is a prominent 
feature of all three "full" , long-wavelength, and pertur- 
bative results, is directly related to the flat background of 
states below the quasiparticle peak in the spectral func- 
tion. The localization-peak feature of the spectral func- 
tion is responsible for the peak in N(lo) at low u>. 




0.04 0.06 

T/2J 

FIG. 9. C M {T) v.s. T for the spin- 1/2 system for x = 
(dashed curve), x — 0.1, and x = 0.2 (dotted and solid 
curves). Dotted curves are the long- wavelength results, solid 
curves are the result of numerical integration using the full 
Green's function (|l0|). The dashed sector shows the 3D 
crossover temperature region T < \fJJi_ for J± = 10 -4 J. 

The calculation of magnetic contribution to the spe- 
cific heat using results for N(ui) and Eq. ( f45[) is straight- 
forward. For the pure system Cm(T) ^T 2 because of 
the two-dimensionality The anomalous density of states 
results into a quasi-linear correction to it. Using the long- 



SC M (T) 



A(x) 



T 



x ln 2 |T/w |+7r 2 /4 ' 



(70) 



where A{x) is a weak function of x. At T > luq (T is 
also in units of Qq) this gives 



5Cm (T) ~ xT ■ const 



(71) 



Fig. M shows our results for the magnetic specific heat 
of the spin-1/2 system (17q = 2 J) v.s. T for x = 
(dashed curve), x — 0.1, and x — 0.2. Dotted curves 
are the results from the long-wavelength expression for 
N(uo), solid curves are from numerical integration using 
Eqs. (EOT) and (Eq). One can see that the results are 
very close and point to the same behavior. The dashed 
sector shows the temperature region T < \fJ~JZ where 
the crossover to the 3D behavior (which provides higher 
powers of T to Cm) should occur. We use the value 
J± = 10 -4 J characteristic for the cuprates. 

Realistically, this picture should be overlapped with 
the phonon contribution to the specific heat. One would 
expect phonons to remain essentially three-dimensional 
even in the layered materials with the characteristic T 3 
contribution to the specific heat at low temperatures 
and thus be negligible in comparison with T 2 and xT 
terms. However, in the case of cuprates the phonon De- 
bye energy is of the order of 400A" pa] which is signifi- 
cantly lower than 2 J ~ 3000A. This makes the phonon 
part of C(T) to deviate from the T 3 behavior already 
at about 20AT, that is around the 3D crossover temper- 
ature for magnetic subsystem. Because of much lower 
Debye energy the specific heat in cuprates is dominated 
by the phonon part |5q ]. Therefore, in order to observe 
the anomalous quasi-linear contribution of the localized 
states to C{T) one needs to use the "reference" mate- 
rial x = and subtract C x= q(T) from the results for 
the systems with x > (we assume that impurities do 
not introduce dramatic changes in the low-energy phonon 
spectra). Another route is to find a quasi-2D system with 
much lower value of J (of the order or less than Debye en- 
ergy for phonons) which would allow a direct observation 
of xT-anomaly from localized states. 

The finite value of the inter-plane coupling together 
with the small anisotropy gaps leads to the finite value 
of the ordering temperature Tjv whose dependence on 
impurity concentration is considered in the next Session. 
The effect of the 3D coupling in the dynamic proper- 
ties, briefly mentioned above in the context of the specific 
heat, is the following. The energy scale introduced by the 
inter-plane coupling r = J±/2J is oj^d = vTr as seen 
from Eqs. (ra), (113) and therefore is rather small for the 
realistic systems of interest (for LCO ui^d ~ 0.01 in the 
units of the magnon bandwidth) . We show in Appendix 
H that the 3D corrections to the 2D scattering are given 
by O(rlnr) which is truly negligible (~ 10~ 4 for LCO). 
Therefore, the only appreciable correction to the dynamic 
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properties from 3D coupling is the low-energy cut-off of 
the logarithmic terms in the self-energy at u> = u^d- As 
we describe in Appendix |] 



P3d{u) ~ -In 

7T 
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at cj -c u 3D , 



(72) 



that is, below the 3D energy scale the real part is a con- 
stant and imaginary part has an extra power of us in 
comparison with the pure 2D form of p(u>). p{ui) remains 
essentially two-dimensional at u> > oj^d- Evidently, this 
proves that the 3D coupling has little or no effect on the 
properties of the spectral functions, dynamical structure 
factor, or density of states at ui > UJ3D- 

However, the 3D coupling does affect some of the lo- 
calization features in the following way . Below the 3D 
energy scale the "stretching factor" ( p5J ) saturates at the 
value a(w3u) and the imaginary part of the self-energy 
acquires an extra power of w. In other words, it should be 
understood as the competition of disorder-induced and 
3D energy scales. Therefore, there are two regions of x. 
First, when x is small enough < x < x* ~ 1/lnr -1 
so that a(w3£>) > 0. In this region the well-defined 
spin waves can be found deep in the low-k, low-w re- 
gion (k,uj <C w 3 d), similar to the quasi- ID problem E7| . 
Concentration x* is defined from the equality of the en- 
ergy scales e~ 7T / 4x = y/r which gives x* ~ 1/ lnr -1 . The 
localization peak in the spectral function at lu ~ lu ~ 
j e -ir/4x j n -j_ ne i ow _ a7; k ^ w region will be replaced by 



\ ' (w) W - 



1 



7r w k a(w 3D y 



at w« uj 3D , 



(73) 



which smoothly vanishes as to goes to zero instead of 
showing a peak. However, the nonlinearity of the spec- 
trum, abnormal damping of the quasiparticles, and the 
flat background of the localized states below cDk are all in 
the 2D-region of k — u> space (u> > lu^d) and will remain 
intact. 

Second region is x > x* where cl(u>3d) < 0- I n this re- 
gion the pole at low-k and low-w becomes pure imaginary 
as in 2D case and the localization peak for low-w, k>u 
reappear above the 3D scale. Above the concentration x* 
all the low-energy excitations are incoherent because the 
2D disorder-induced energy scale ujq (localization length 
£) is larger (shorter) than the 3D energy scale w^d (length 
scale 1/\/t) so the spin waves lose their coherence before 
they can propagate in 3D. A self-consistent calculation is 
required to determine accurately the value of x* and the 
details of the 3D to 2D crossover. Our estimation gives 
x* ~ 0.1 - 0.2 for r - 10~ 4 . 

Thus, we find that the 3D coupling for the realistic 
materials will modify the 2D density of states, structure 
factor, and specific heat only at the energies (tempera- 
tures) u < wsd — 0.01 and at impurity concentrations 
x < x* ~ 0.1 — 0.2. The estimated value of the 3D cou- 
pling r c which would make x* larger than the percolation 
threshold is r r ~ 0.01. 



The consideration given above also applies to the case 
of small anisotropics introducing gaps in the spectrum 
with a modified r = r e ff accumulating the total effect of 
the gaps and 3D coupling. It should be noted that the 
incoherence comes from the averaging procedure which 
converts the dissipation of momentum into the dissipa- 
tion of the energy. Therefore, the overdamped excitations 
should be understood as diffusive. It is interesting that 
it requires 2D and "strong" disorder to restrict the num- 
ber of Euclidean paths for spin waves and to break down 
the description of the problem in terms of an effective 
medium. 



IV. STATIC PROPERTIES 

The static properties such as average staggered magne- 
tization M(x,T), Neel temperature T/v(x), and 2D cor- 
relation length £(T, x) are considered in this Section. 

The average on-site magnetic moment Eq. (46) for 
randomly diluted AF with the averaging over magnetic 
sites M(x) = ^i\Sf\/N m , see Q, can be exp ressed 
through the integral of the spectral functions ( |4l|) as: 



M(x, T) = S-A- SM(x, T) , 
5M(x,T) = ^ r*B{») du. 



(74) 



k 



Wk 



[< k (a;)-7k<, k (o 



where A = X^k^k — 0.1966 is the zero-point spin devia- 
tion, ris(cj) = [e u / T — 1] _1 is the Bose distribution func- 
tion, subscript R denotes retarded. Note that one should 
not expect this formula to be valid at large doping level, 
x close to x p , since our approach neglects decoupled clus- 
ters and interactions of impurities. However, at not too 
large x these effects should be negligible and one expects 
Eq. ([74]) to be adequate. We would also like to note 
here that our definition of M(x,T) is physically equiva- 
lent to the "quantum-mechanical factor" of the averaged 
staggered magnetization, the definition used in the re- 
cent Monte Carlo study pq] . In other words, the "classi- 
cal" ("geometrical") effect of dilution on magnetization, 
which simply accounts for the decrease of the magnetic 
substance, is multiplicative to the quantum effects and is 
not taken into account in Eq. (u4T). 

First we address the question of the presence of explicit 
divergences in the integral Eq. ( |74| ) which would point to 
the instability of the long-range order discussed in Refs. 

@0- At T = ° ns H = -0(-<*>) and the impurity- 
induced quantum reduction of the magnetization, which 
can be interpreted as a result of the "condensation of 
magnons" , is given by 



6M(x) 



k 



7-1 w k 



M] , (75) 



where we use that the spectral functions are zero out- 
side of the magnon band lu 2 > 1 . Since the perturbative 
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result Eq. (|5l]) suggests the instability at small wave- 
vectors the long-wavelength expressions for the spectral 
functions can be used for our analysis. From the form of 
the spectral functions in Eqs. (|5J), (|6l| ) one can readily 
see that the integral over u) is always finite. The integra- 
tion over k is two-dimensional but has a factor of 1/wk 
in the integrand. From our expression of the spectral 
functions in the intermediate and localization peak en- 
ergy ranges Eqs. (|58|), (p9|), ( J62] ) one may suggest that 
there is another 1/c^k in the integrand which would lead 
to the logarithmic divergency. However, these expression 
are obtained by neglecting uj 2 in comparison with w k and 
thus are valid at w^ ^$> lu only. At lower k the convergence 
of the integral is restored. To show that more explicitly 
one can use the ^-expanded form of the Green's function 
Eq. ( |67| ) for A^^lu) and an equivalent expression for 

4? k M 



Carlo data from Ref. |2q| (filled circles), and NQR data 
(open circles) from Ref. |36| are also shown. Note that 
the original Monte Carlo data of Ref. p5| are normalized 
by the total number of sites while both NQR and our 
results are averaged over the magnetic sites only. In or- 
der to extract the same quantity from the Monte Carlo 
data we divided them by the classical probability to find 
a spin-occupied site within the infinite cluster |48|| . A re- 
cent Monte Carlo study Ref. (26) provided an analytical 
expression for the fit of the "quantum- mechanical factor" 
in the magnetization (see the comment after Eq. (J74|)) 
which we plot in Fig. [l(J as well (dashed line). One can 
see a very good agreement of our results with numeri- 
cal data up to high concentrations. The oxidation of the 
crystals can be the reason of a faster decrease of M(x) 
in NQR data. 



G k >)^Gr>)£MGT>) 



(76) 



Since all E's are linear in x this provides an expression 
for the linear in x term in the staggered magnetization: 



SM(x) ~ xB 



E 



° _dw 

-1 7TCo> k 



(u) - LUk) 2 uj 2 - w k 



+ 



7k Re£ k 2 ( Wk ) 



k 



2ul 



(77) 



In the long-wavelength limit this gives 



8M(x) ~xB~-Y I c 
n k ^o 

X sr^ 1 

+ -> — 



1 



(lu + lu^) 2 UJ 2 



In 



L^k 

4 



-7T 2 /4 



(78) 



where the strongest divergency of the integrand is In k dk 
and all integrals are convergent. 

Numerical integration of the expression in Eq. (|77| ) 
without the long-wavelength approximation gives the 
suppression rate of the staggered magnetization M(x) ~ 
M(0) - Bx with B = 0.209(8). For S = 1/2 it 
gives the slope of the normalized staggered magnetiza- 
tion M(x)/M(0) ~ 1 - Bx/(S - A) ~ 1 - 0.691(5) • x. 
It is interesting to note that the second Born approxima- 
tion to the impurity scattering gives three times smaller 
rate Bbotu — 0.0725 showing the necessity of the full 
T-matrix treatment of the problem. The estimation of 
B, given in the previous study Ref. [ p7[ using 1/z ap- 
proximation for the expression similar to our Eq. (77), 
provides even smaller B 1 / z ~ 1/z 2 = 0.0625 showing yet 
another inadequacy of that work. 

We have also performed a numerical integration in Eq. 
( |75| ) for the impurity-induced reduction of the staggered 
magnetization without x-expansion. This yields the re- 
sults presented in Fig. nfl for S = 1/2 (solid line). Monte 
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FIG. 10. Average staggered magnetization v.s. x. Our 
results from Eq. ( pq ) (solid line), Monte Carlo data (open 
circles, Ref. ||), NQR data (filled circles, Ref. @), and the 
fit of Monte Carlo data from Ref. fE6l are shown. 



The absolute value of impurity-induced quantum fluc- 
tuations 8M(x) is independent of S in the linear spin- 
wave approximation similar to the quantum reduction 
of S by zero-point fluctuations A. We plot our results 
for 5M(x) in Fig. [I]] in order to emphasize the agree- 
ment with the MC data for S = 1/2 (circles) and 5 = 1 
(squares), which show only weak S'-dependence. 

It is worth mentioning here that the discrete static 
quantities, zero-point spin deviations at the neighboring 
sites around impurities in an AF, were studied using spin- 
wave theory and Green's functions methods since sixties 
|57j with most recent results obtained in Refs. 0, [p3[ . 
Quite remarkably, these results |5q| were found to be in a 
very good agreement with the recent Monte Carlo studies 
of 2D S = 1/2 Heisenberg model with impurities, Ref. 
||. Note that while Refs. | p3| , J57J] were focused on the 
discrete quantities our results concern the averaged ones. 



15 



0.25 

0.20 

■0.15 



o MC, S=l/2 
■ MC, S=l 
theory 




FIG. 11. The absolute value of SM(x) from Eq. @ (line) 
and Monte Carlo data, Ref. [|25[, for S — 1/2 (circles) and 
5 = 1 (squares). 

At T > Eq. (|7J) for the staggered magnetic moment 
can be rewritten separating the quantum, T = 0, and 
thermal, T-dependent, parts 



M(x, T) = S - A - 5M{x) - SM 1 (x, T) , 
^ Jo ^k 



(79) 



M 



-2 7k A k 2 M] , 

where 8M(x) is the zero-temperature part given in Eq. 
( |75| ) and we used evident symmetries of the spectral 
functions with respect to u) — ► —uj and that %(u) = 
-1 - ub{— w)- 

For the true 2D system at x = and T > thermal 
fluctuation destroy the LRO which manifests itself as a 
log-divergency of the thermal correction to the magneti- 
zation 



SM T {0,T) 



E 

k 



ub(u)) dco 

w k 

T Tduj 



S(UJ - LU k ) 



(80) 



where we use that ns(u) ~ T/lj at T <C w. The 3D 
coupling provides a cut-off to this divergency in a quasi- 
2D problem which yields the finite value of the thermal 
correction 



SM T (0,T)~-[ n B (w)dcj~ -Tin 



(81) 



and the finite value of the Neel temperature whose mean- 
field value can be found from the condition M(0, T) = 

= S - A - 8M T (0, T) which gives 



-N 



tt(S-A) 
lnr- 1 



«1, 



(82) 



in units of 45 J. T/v vanishes when r — » 0. 

One would expect that the thermal part of the stag- 
gered magnetization for the diluted system may possess 
other divergences, stronger than the simple log-w for the 
pure system. In fact, this suggestion is quite natural since 
the spectrum is not linear and, therefore, the non-linear 
corrections must show themselves up. Indeed, since the 
correction to the spectrum is <5w k ~ xw k ln|w| (51) one 
immediately suggests that the thermal part of the mag- 
netization should acquire a term 



InM^ 



(83) 



However, we show that such anomalous terms from diago- 
nal and off-diagonal spectral functions cancel each other. 
As a result, there is no signature of any new divergency 
in this quantity caused by the anomalies of the spectrum. 
Using the x-expanded form for the Green's functions 
(pTJ) , (um in the long- wavelength approximation one finds 
the diagonal 



^J(x,T) =£-L(a k a k )^W 



Ub{ijj) diu 



X < 7T<5(w — Wk) 
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parts of the temperature dependent 8M T (x,T), where 
we kept only C(ln \uj\) and 0(1) terms in the integrand, 
p'{uS) = Rep(w), integration by parts was used in 6Mj, 
superscript T in the averages means the thermal part. 
The total result is 



SM T (x,T) 



1 + X 7T 



SM T (0,T) 



(86) 



which shows that the thermal correction is enhanced 
by impurities but there is no new divergency associated 
with them in this quantity. Suppression rate of the Neel 
temperature can be readily obtained from the condition 
M±x,T\ = = 5 - A - 5M(x) - SM(x,T) using Eq. 
( |S6| ) which gives: 
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For 5 = 1/2 this gives A 1/2 = 3.196(5) and for 5 = 5/2 
it is A 5 / 2 = 2.600(4). It is important to note that these 
suppression rates point to x c (l/2) ~ 0.31 and x c (5/2) ~ 
0.38, both below x p , so that one may suggest that in order 
to have the phase transition at the classical percolation 
threshold the Tjv(x)/Tjv(0) curves should have a rather 
unusual concave form. 

It is interesting to compare our result for the decline 
rate of T/v (x) ( p7| ) to the answers of different approaches 
to the same problem and to the results for similar mod- 
els. A naive mean-field treatment of the impurity ef- 
fects as simple renormalization of magnetic coupling gives 
Tjv(x)/Tjv(0) = 1 — x. Application of our formalism to 
the Ising limit of the 2D problem gives Tv(x)/Tv(0) = 
1 — A 1 x with A 1 ~ 1.37 (see Appendix |FJ) which is very 
close to the RPA answer Aj^p A = 1.33 and below the ex- 
act answer A\ xact ~ 1.57 J58|. For the 2D Ising magnets 
Tjv(x) v.s. x has a more traditional convex form p7[ . 
Previous result on the suppression rate of Tv for the 2D 
Heisenberg model [|9| is T N (x)/T N (0) = 1 — ttx which is 
obtained using Green's function technique and spin- wave 
theory in approximations very similar to ours. However, 
Ref. |p9| misses ~2/n and neglects 1/5 terms. 
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FIG. 12. Tn(x)/T n (0) v.s. x for S = 1/2. Results of 
numerical integration in Eq. (174) (dashed line), analytical 
linear- a; slope (1 — Ai/ 2 %) Eq. (87|) (solid line), /jSR (dia- 
monds) pij and magnetic susceptibility (circles) jl9| data for 
Zn-doped LCO, and ESR (squares) |22J of Zn-doped copper 
formate tetrahydrate Cui-^Z^Mg^HCOaJ-HaO. 

We have also performed a numerical integration in Eq. 
(JM) and solved an implicit equation M(x, T/v) = on 
Tn(x) numerically. This procedure requires the finite 
3D coupling and the use of quasi-2D form of the spectral 
functions. Since the integration involves an additional 
dimension and the 3D region is quite narrow the conver- 
gence of the result as a function of number of k, w-points 
at small x can be an issue. We plot our numerical results 
for T N (x) /T N (0) for the case of 5 = 1/2 in Fig. || to- 
gether with the analytical slope Eq. ( p7| ) with Ax/2 = 3.2 
and experimental data. Experimental data are obtained 



by /iSR plfl and magnetic susceptibility measurements 
|l9| of LCO systems and by ESR [22J of Zn-doped copper 
formate tetrahydrate, a layered quasi-2D AF. One can see 
that our linear-x results agree very well with the exper- 
imental data up to a rather high doping level x ~ 0.25. 
There is a slight disagreement between our own numerical 
and linear- a; analytical results already at small x which 
may be connected not only to the numerical accuracy 
but to the corrections of the order ~ xT^/J ~ xj In t -1 . 
Note that the linear- a: result is free from such corrections 
since it is obtained in the r — ► limit. 

As it is discussed extensively in Ref. ]6(J the spin- 
wave theory for layered materials is not really adequate 
at T ~ Tv because of the lack of the kinematic con- 
straints. When it is applied to the mean-field equation 
M(x, T/v) = it tends to overestimate the absolute value 
of the Neel temperature and has some other artifacts such 
as M(T) - T N - T at T - T N . This may also provide 
an additional x-dependence in our numerical values of 
T N {x)/T N (0). 
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FIG. 13. T N (x)/T N (0) v.s. x for S = 5/2. Analyti- 
cal linear- a; results (1 — A 5 / 2 x) Eq. ( p7[ ) (solid line), and 
(circles) J24J magnetic susceptibility, specific heat data for 
Mni_ I Zn(Mg,Cd),(HC0 2 ) 2 2(NH 2 ) 2 CO. 



shows our analytical slope for Tn(x)/Tn(0) 




Eq. 1871) with A 5 / 2 = 2.6 for the case of 

5 = 5/2 and experimental data from mag- 
netic susceptibility and specific heat measurements of 
Mm_ a; Zn(Mg,Cd) a; (HC02)22(NH2)2CO, a layered 5 = 
5/2 material, [E4J. One can see that while the scatter- 
ing of experimental points seem to be smaller than in 
5 = 1/2 case the linear- a; result fits them very closely up 
to x = 0.2. The older T N (x)/T N (0) data for the more 
traditional 5 = 5/2 material K 2 Mni_ :E Mg :E F4 J23| show 
a big scattering of the data which allows almost any rea- 
sonable fit @. 

It is worth mentioning that the numerical results of 
our approach for T/v(x) bends inward at larger values of 
x and show the above mentioned concave form, which 



17 



has been recently observed experimentally for LCO com- 
pounds |L8| and has been anticipated in other works J24| . 
While our approach is certainly not adequate at such 
high impurity concentrations and tends to overestimate 
the value of Tjv(x) in comparison with experiments, it 
nevertheless points to the same physics. We interpret 
this behavior as due to localization effects which tend to 
reduce the role of quantum and thermal fluctuations in 
the destruction of the long-range order. 

In the paramagnetic phase above the Neel ordering 
temperature the 3D coupling is irrelevant and the spin 
fluctuations in the layered AF system are characterized 
by the in-plane correlation length £20 which is exponen- 
tially diverging with 1/T as T — > 0. The correlation 
length is uniquely determined by the T — properties of 
the system such as spin stiffness constant p s . 

The correlation length can be derived from the mod- 
ified spin- wave theory, as was suggested in Ref. pQ ], by 
introducing a chemical potential for magnons which pro- 
duces a gap in the spin-wave dispersion and then by re- 
solving a constraint (Sf ) = which defines the correla- 
tion length sclf-consistcntly. The result of such calcula- 
tions at x = is fbOl 



£(T)^-exp( — 



(88) 



One should bear in mind, however, that while this 
approach gives the correct exponential behavior of 
£,2d(T) it provides a prefactor equivalent to the one- 
loop renormalization-group result [|l|.[2|. This prefac- 
tor must be modified according to the higher order 
renormalization-group treatment [pl| which gives |62| 



ecn 



exp 



f2irp 6 



2(4^p s + T) ^\T 



(89) 



This expression shows excellent agreement with experi- 
ments and Monte Carlo data j|j2f|,^j . This discrepancy 
between the results of modified spin- wave theory and re- 
sult of more exact, non-perturbative approach is of the 
same origin as the overestimation of the Tn by the mean- 
field solution of (Sf) — equation |p0| . 

We generalize the approach of Ref. j50| for the case of 
an AF with impurities and obtain for the constraint: 



S 



2 ^ 

k 



Wkfa) 



(90) 



4- Wkfr)) i 



«|c«k) 



7k(44) 



where 0J\^{rj) = y\L — ?7 2 7^ and magnon averages are 
given by the integrals of the spectral functions A^(uj), 
A^{uj) from Eqs. (J40|), (fllj) in which the gapped form of 
the spin- wave spectrum is used. 

We have performed a numerical integration in Eq. 
( pp| ) and calculated the correlation length S,2d(x,T) — 
rj 2 / \J&{\ — r] 2 ) as a function of T for several values of 



x. We fit the results of such a numerical procedure 
in a wide temperature range almost exactly with the 
help of original Takahashi formula, Eq. ( p8[ ) , with spin- 
stiffness ps(x) being a free parameter. These fitting val- 
ues of p s (x)/p s (0) v.s. x follow closely our result for 
Tjv(x)/Tjv(0) dependence, Fig. [l^. Such a result can 
be anticipated from the mean-field picture of the or- 
dering in layered systems. The transition occurs when 
the inter-plane coupling is strong enough to stabilize 
the LRO in comparison with the thermal fluctuations: 
Jj_M 2 (x)£_ 2 (x,T N (xj) « T N . If the correlation length 
preserves its exponential form the dominant part of the 
left-hand side comes from e 2 ^P=>{x)/T N (x) anc j one j mmc _ 

diately arrives to 



Ps(x) = 

P.(0) T N (0) 



Tjv(x) 



+ 0(x/lllT- 1 ) . 



(91) 



Therefore, the important conclusion one can make from 
our analysis is that (i) the correlation length should fol- 
low the x — type of behavior Eq. (|8|) with x-dependent 
p s , at least for not too low T and not too high x, (ii) 
p s (x)/ Ps (0)-T N (x)/T N (0). 

Our Fig. [14| shows a semi-log plot of £(x, T) given 
by formula in Eq. (|89| ) with p s (x) — p s (0)(l — Ai/ 2 x), 
A1/2 is from Eq. (|87|), v.s. J/T for x = (dashed line), 
x = 0.1, x = 0.2, and x = 0.3 (solid lines). An impor- 
tant observation can be made here. At small x 2irp s is 
of the order of J and at all reasonable temperatures the 
dominant behavior is exponential in J/T (straight line 
in the semi- log scale). When the spin stiffness becomes 
small (p s <ti J) there is an additional temperature range 
J 3> T 3> p a where the exponential behavior is not seen 
yet while the prefactor gives a log( J/T) behavior of the 
log(£) clearly seen for x = 0.3. The experimentally ob- 
served deviation from the simple exponential behavior of 
the correlation length £(T, x) v.s. 1/T |Q can be related 
to this effect. 




FIG. 14. £(x,T) from Eq. (|8S|) 

with p a (x) = p a (0)(l - A 1/2 x), from Eq. @, v.s. J/T 
for x — (dashed line), x — 0.1, x — 0.2, and x = 0.3 (solid 
lines) . 
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At larger doping level close to the percolation thresh- 
old one expects a new length-scale to appear. This 
length-scale is associated with the crossover from trans- 
lational invariance to the self-similarity in the percola- 
tive systems |64| . Below x p the "geometrical length" , 



£g oc I a; 



~ u , separates the regions of Euclidean 



and fractal geometry. Above x p , where no infinite clus- 
ters left, £g is the characteristic size of the finite clus- 
ters. Earlier experimental studies of the 2D and 3D Ising 
(Rb 2 Co 2; Mgi_ 2 ;F4 and Fei-zZna^) [p5p^ ] and near- 
Heisenberg (Rb 2 Mn £C Mg 1 _ a .F 4 and Mn^^Zn^Fj) [^JG8| 
systems close to x p have demonstrated that the static 
structure factor <S(q) contains contributions from both 
"thermal" and "geometrical" lengths in agreement with 
the theoretical studies [|3 2||6 Sf| . The experimental data 
suggest that these lengths combine in the simplest possi- 
ble form £ _1 = £ _1 (T) + S.^ix) and that the Lorentzian 
form of the structure factor near ordering vector is pre- 
served. Yet another interesting result of the proximity 
to the percolation is that at x < x p below T/y the elas- 
tic Bragg peak should be accompanied by the Lorentzian 
whose width at T = is solely defined by the inverse " ge- 
ometrical length" £q . One would expect similar effects 
to be observed in newly available LCO systems close to 
;c p |Il. 

It is not clear, however, whether the localization ef- 
fects in the infinite cluster, which we discuss in this work, 
can manifest themselves in the static structure factor 
or correlation length. Such contributions, if they exist, 
may lead to a new behavior of correlation length, dif- 
ferent from the simple renormalization of spin stiffness. 
However, in our approach the potential sources of such 
anomalous terms appear in the higher order in x (~ x 2 ) 
and, most certainly, do not affect the results for the ex- 
perimentally reachable domain of lengths £ < 200a above 
which the ordering occurs. At larger concentrations x 
such contributions can become important for shorter cor- 
relation lengths but in reality they might be screened by 
the similar effects from the decoupled clusters. 

Theoretically, it is very intriguing if such localization 
effects of the infinite cluster can really affect the behavior 
of correlation length. We reserve this subject for the 
further study. 



V. CONCLUSIONS 

We have studied the problem of diluted 2D and quasi- 
2D quantum Heisenberg antiferromagnets in a tetrago- 
nal lattice making use of linear spin-wave theory and T- 
matrix approach. We have shown, contrary to the earlier 
findings, that the 2D is not the lower critical dimension 
for this kind of disorder and that at T = long-range 
order persists up to concentrations close to the classical 
percolation threshold. These results are consistent with 
Monte Carlo simulations in large lattices p5[ . In agree- 
ment with earlier works on this subject, which studied 



the problem in the leading order of the dilution fraction x 
p0| , p7| , we found that the spin- wave spectrum is strongly 
modified by disorder. However, contrary to these works 
we have shown that this result does not imply an insta- 
bility of the system to a paramagnetic phase. It rather 
indicates magnon localization on a length scale £, expo- 
nentially large in 1/x. We have shown that this new 
length-scale appears explicitly in the dynamic properties 
such as the dynamical structure factor <S(k, to) , Eqs. (p7|)- 
(pq), which can be measured directly in neutron scatter- 
ing experiments, and the magnon density of states N(u), 
Eqs. (pq), ( p9[ ) which is directly related to the mag- 
netic specific heat Eqs. (|70|), (frlj). The measurement of 
such quantities will provide a direct test of our theory. 
Furthermore, we show that the static properties such as 
the zero-temperature staggered magnetization M(x), Eq. 
( |74| ) and Neel temperature Tjy(x), (in the quasi-2D case) 
do not show any anomaly associated with the spectrum 
and are finite up to the concentration close to the classi- 
cal percolation threshold. These results are in a quanti- 
tative agreement with the NQR || 5 /iSR @, ESR ||, 
and magnetic susceptibility |L9] measurements in differ- 
ent compounds as well as with the Monte Carlo data [B5| . 

We have shown that the effect of dilution of an AF with 
non-magnetic impurities is quite strong because dilution 
removes completely spin degrees of freedom from the im- 
purity site and, therefore, the spin waves are strongly 
scattered. Moreover, the low dimensionality of the sys- 
tem constraints significantly the phase space for scatter- 
ing leading to the localization effects. We have shown 
that the hydrodynamic description of the problem breaks 
down for length-scales larger than £ and the spin excita- 
tions become diffusive instead of ballistic. The conven- 
tional averaging procedure which is used to treat disorder 
does not lead to an effective medium with renormalized 
parameters. Therefore, one needs to use a different ap- 
proach for length-scales larger than £, the problem which 
is beyond the scope of this paper. 

In fact, the physics of localization described in our 
work has similarities to the Anderson localization for 
non-interacting electrons in disordered lattices where the 
statistics of the excitations does not matter JT0| . Note 
that our problem should be close to the problem of lo- 
calization of relativistic bosons (with chemical potential 
\x = 0) in a random potential. On the other hand, that 
problem is related to the problem of disorder in Bose- 
Hubbard model where non-relativistic bosons with ki- 
netic energy J interact through the local Coulomb term 
U J71J . In the latter model the Bose glass phase appears 
for small J at zero chemical potential, and transition into 
a superfluid state is possible when J is large enough. In 
our case superfluidity is not possible but we may con- 
jecture that our localized phase is somewhat similar to 
the Bose glass phase and magnons are trapped in the re- 
gions which are more ordered than in average. It is not 
clear, however, if the relativistic nature of the bosons is 
important for the nature of localization. 

Furthermore, we find the close similarity of our prob- 
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lem to the problem of disorder in 2D d-wave supercon- 
ductors |44|,[52|] ■ The large enhancement of the density of 
states at low frequencies in our case, which comes about 
because of the redistribution of spectral weight over the 
entire Brillouin zone, is reminiscent of that problem. In 
rf-wave superconductors the elementary excitations are 
nodal quasiparticles, or relativistic (Dirac) fermions. It 
is known that for these excitations localization occurs on 
a length scale £l (localization length) which is an expo- 
nential function of the conductance a: £l oc e a l a ° |72| ] 
where er — e 2 /h. Since in the dilute limit one expects 
the conductance to diverge with x (that is, a ex 1/x) the 
localization length has the same type of non-analytic de- 
pendence on x as in our case. However, it is not clear how 
(if possible) the two problems map onto each other. A 
further investigation, beyond the scope of this paper, can 
clarify the connection of the diluted antiferromagnet with 
other similar problems of disorder in low-dimensional sys- 
tems. 

In summary, we have presented a comprehensive study 
of diluted quantum Heisenberg antiferromagnets in 2D 
and quasi-2D. We have shown that while the dynamic 
properties possess anomalies associated with magnon lo- 
calization the static properties are free from such anoma- 
lies. Thus, in low-dimensional systems with disorder the 
connection between static and dynamic quantities is not 
straightforward. We have compared our results to the nu- 
merical simulations and experimental data with a very 
good agreement. We have also proposed other experi- 
ments which can further test the results of our theory. 
Altogether this provides a self-consistent picture of the 
effects of disorder in low-dimensional quantum antiferro- 
magnets. 



Vk!,k 2 = I dri dr 2 (j}i 1 (ri)V rur2 (j) k2 (r2) (Al) 

dri / dr 2 ^>k 1 (ri)C / /K 1 ,r 2 ^0k 2 (r2), 



where Ui is any symmetric operator in the group of 
tetragonal symmetry, and k (r) is a plane wave func- 
tion, 0k (r) = (2tt) ' e lkr , which can be decomposed by 
projection operators: 



p 71 — 1 



where 



0k (r) 



(*>) 



9 t-? 



\toh 



(A2) 



(A3) 



where the set of functions, {0k(r)n }n=i i , form a ba- 
sis of the p th irreducible representation, and l p is the 
dimension of p th irreducible representation; D^ p \Ui) nn 
is the diagonal matrix elements of the p th irreducible 
representation for the symmetric operator Ui in point 
group Dih whose order is g (= 16). We readily project 
the potential into irreducible representations as Vk ll k 2 = 

E K (p 

£-^p ki 



(p) 



k 2 ' 



where 



n 






dk 3 / dk 4 A ki k3 Vk3 :k4 A k4jk2 , 
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APPENDIX A: TETRAGONAL LATTICE GROUP 
THEORY 



where 



We first resolve the scattering potential V 
r-space by inserting the closure relations 173] 



Ofiin 



^kk 2 = / dr0 ki (r)f/^ k2 (r). (A5) 



Using the tetragonal symmetry group one notices that 
each A l ki k2 is a <5-function. Thus, the scattering po- 
tential V^J k2 Eq. (|9|) can be decomposed into channels 
of irreducible representations. The non-zero orthogonal 
channels are (before the Bogolyubov transformation): 
Aip (s-wave): 



K 



|ski; 



'ki,k 2 
E M (in-plane p- waves 



> s k 2 



"ki,k 2 

Big (d-wave): 



Pki 



(y)\ 



'ki/ 



(Pk^ 



V, 



ki,k 2 



|d kl )«)(d k2 | 



(A6) 



(A7) 



(A8) 



A 2 « (p z -wave): 
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V, 



ki,k 2 



|Pki) ® (P: 



k2 I ' 



(A9) 



where (s k | = [1, 7k ], <s^| = ^[l,^], <Pk to) l = 
[0, 1] sinfe x(3/) /v/2, (d k | = [0,1] 7^ and (p k | = 
•Jr sin fe z [0, 11 /y2. Bogoliubov transformation yields 
Eqns. @-©. 



APPENDIX B: 3D T-MATRIX 

In this Appendix we provide the solution of the s- 
wave T-matrix equation in the tetragonal lattice for the 
arbitrary relative value of the inter-plane and in-plane 
exchange integrals t = J±/2J. With this solution we 
demonstrate the smallness of the 3D corrections to the 
2D result in the quasi-2D case (r <C 1). 

After some algebra one can solve the T-matrix equa- 
tion ( p5| ) with the s-wave scattering potential from Eg . 
( |l8| ) (sublattice A), Uk, Uk, and Wk from Eqs. p2), pa ) 
and obtain: 



r^M = |ak>®<«k'| -TxH 



+ (k> 



r 2 (w) 



(Bl) 



(41 



(sk'|) -TaH, 



where the w-dependence of the in-plane scattering (first 
term) is given by expression which is formally similar to 
the pure 2D result Eq. ©: 



TiH 



(1 + u))p(w)+tRi(w) 



1 



lu 1 — lu(1 + lu)p(lu) + tD(u>) 



(B2) 



with p(u>) given by Eq. (E9). Note that the integration 
over p in this case is three-dimensional. The inter-plane 
scattering is: 



r 2 (w) = -t + — 



t 2 R 2 (lu) 



lu 1 — u(l + uj)p(uj) + tD(u>) 
The w-dependence of the cross-term is given by: 

tR 3 {lu) 



T 3 (lu) = 



lu 1 — w(l + uj)p(u>) + tD(lu) 



(B3) 



(B4) 



All three parts of the scattering matrix possess the same 
"unphysical" 1/lu contribution discussed in the text. Ap- 
plication of the projection procedure Eqs. (|30|), ( [H| ) to 
this problem is out of the scope of this Appendix. 

The auxiliary functions D and Ri are given by rather 
cumbersome combinations of u>, p(u>) and two additional 
integrals: 



*M = E 



(7^) 2 ,. A v- 7 P 7p~ 



£M = E' 



(B5) 



with 7 P , 7^ , w p from Eqs. (|), @. Note that at r -> 
a(w) -> p(w)/2 and /?(w) -> 0. 



The expressions for D and i?i are: 
D = P - LuP 2 , 

fli = P - a + P 2 , (B6) 

R 2 = a — (w - t)P 2 , 
i? 3 = (P + a - P)/2 + r(p - o)/2 + P 2 , 

where the following shorthand notations are used 

P = %(p + a)-2(3, (B7) 

Pi = 7o P a — P + a — wP — u) pa . 

There is no assumption on the value of r made in these 
formulae. 

At t <C 1 and w«l(u can be still » r) one can show 
that: 



Dc 

R-2 



3p/2, 
-p/2, 



i?3 



-P 2 /2 + P 
0(rp 2 ) . 



(B8) 



In the same limit r < 1 and w <C 1 the w-dependent 
parts of the scattering matrix become (we simply omit 
the unphysical 1/lo terms here): 

,-r(p 2 



Ti(w) ~p-r{p A - p) 

T 2 (uj)~-T + T 2 p/2 , 

r 3 (w) = o(rV) ■ 



(B9) 



Recall that in 2D Ti(w) ~ p and T 2 (uj) = T 3 (uj) = 0. 
Since Rep ~ In |w| at w> ^Jt and Rep ~ In \r\ at u < 
\[At the largest relative correction to the 2D terms in the 
scattering matrix is C(rln(r)). The same statement can 
be proved for all higher powers of lu in Eq. ( p32j ) without 
making lu <C 1 assumption. 

The conclusion is, once again, that at t < 1 one can 
safely drop all terms explicitly proportional to r in Eqs. 



(B2)-(B4) and thus arrive to the purely 2D expression 
for the scattering matrix given in Eq. (Eq). The only 
modification in the quasi-2D case versus 2D case is the 
change of the behavior of p(u>) at low lu, whose real part 
saturates at lu < v4r and imaginary part acquires an 
extra power in lu (see Appendix ph . 



APPENDIX C: ELLIPTIC INTEGRALS 

The energy-dependent part of the T-matrix Eqns. 
(Eq), (28) is expressed through the integrals of the 
Green's functions Eqn. (E9). These integrals can be eval- 
uated in the case of 2D and are given by combinations of 
complete elliptic integrals of the first and second kind: 



PM = E -tz 



W - LU* 

p p 



2 

TTLU' 



, x V- (7 P ) 2 
p d (u>) = }_^ 



LU Z — LUp 



K(Lu') + iK(iu) , (CI) 



lu 2 K(lu') - 2E(lu') 



i( (lu 2 -2)K(lu)+E{lu) 
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where u>' — y/1 — ui 2 , K and E are the complete elliptic 
integrals of the first and second kind, respectively |74| . 
In the low-energy limit: 



p{u) = - In |w/4| - $ , 

4 
/Od(w) = 1 . 



(C2) 



APPENDIX D: PROJECTION OF UNPHYSICAL 

STATES 

After introduction of the fictitious magnetic field to 
project out the unphysical on-site mode the s-wave scat- 
tering potential (sublattice A) is given by the sum of two 
terms from Eqs. @, @: 

^A^total = _ |gk) ^ (sk/ | + ^| Ask) g, (Ask ,| ; 

with (s k |=^k[wk, -Uk] , (As k | = [u k , u k ] . (Dl) 

One immediately suggest the form of the solution of the 
T-matrix equation: 

*££(") = l*k> ® <«k'| -TiH 

+ |As k )®(A Sk ,| • r 2 («) (D2) 

+ (|As k ) <g> (s k '| + |s k ) ® (Ask»|) • r 3 (w) , 

and, after some algebra, one finds: 

H z (l + u)p(uj) + l 



TiM = 



r 2 («) = - 
r 3 (w) 



[l-w(l + w)p(u)](H z -w) 

luH z 



H z -w 
H z 



(D3) 



H z -u> ' 



which yield the answer given in Eqs. (p2|)-(|34|) in the 
limit H, — > oo. 



APPENDIX E: 3D p(w) 

The key ingredient of the low-energy T-matrix scatter- 
ing is given by the integral of the Green's function over 
k, p(w) (p9|). Appendix O gives an analytical expression 
of p(uj) in the 2D case. In the quasi-2D case the inter- 
plane coupling provides a cut-off in the logarithm and 
gives an extra power of u) in the imaginary part of the 
integral in the 3D energy range. This can be obtained 
explicitly using 3D form of the spin-wave dispersion Eq. 

(PD^k 



il 



V¥ _ 

In the limit -/F = a/Jj_/2J <C 1, and w < 1 (for 
arbitrary ui/^/t) one obtains for the real part of p{ui): 



Rep(w) = — In 

7T 



Rep(u>) — — In 

7T 



+ 0(r,u/), for lo< V4t , 



+ a/w 2 - 4r 



+ G(r, W 2 ), (El) 

for w > V4r , 



at (j > v4r the 3D energy scale is irrelevant and 
Rep(w) = — In |w/4| is back to its 2D form. Imaginary 
part of p(w) is 



ImpH = -i arccos N^ + r)^-^ -l\ + ^ 

= arccos ( 1 - 77- J +C(r 2 ,w 2 ), 

7r V 2r / 



for w < Vir , (E2) 



Imp(w) = -1 + 0{lo 2 ), for w > V4r , 



at small w <C v 4t deep into the 3D range of energies 

■k-Jt 



Imp(uj) — —-^ is linear in lc. 



APPENDIX F: T N ( X ) FOR THE ISING PROBLEM 

In this Appendix we apply the formalism of our work to 
the problem of Tjv(x) v.s. x dependence for the Ising S — 
1/2 case. While the spin- wave approximation is much less 
adequate in the Ising limit than for the pure Heisenberg 
model it is nevertheless a very instructive exercise which 
gives a quantitatively correct answer. 

Quadratic part of the 2D S = 1/2 Ising model in the 
spin-wave approximation reads as 



U_ 
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Ho + Hi- 



(Fl) 



^4^^^"-*^^, 



i,k,k' 



with 



Vk,k' = 7k-k' 



(F2) 



where we omit from the beginning the "unphysical" term 
which will result in 10 — mode. T-matrix gives the total 
result for all scattering channels: 



tk'M— 7^ 



u- 1 



w - 3/4 ' 
where we used the property 

X^k- P 7p-k' = 7k-k'/ 4 ■ 
p 

The self-energy is then given by: 



(F3) 



(F4) 
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£(u>) = — x ■ 



1 



u - 3/4 ' 
The Green's function has two poles now: 

G[W) -^7~1 a; -3/4 + a; ' 
and the spectral function is given by two <5-peaks: 
1 



(F5) 



(F6) 



A{u) 



■ [5{uj - 1) + 4x5(w - 3/4 + x)] . (F7) 



1 + Ax ' 

Neel temperature is defined from the condition: 



(S'XTjv,*) 



dw n B (w)4(w) = , (F8) 



which transforms to 
l + 4x 



n B (l) + 4a;n B (3/4-a;) 



(F9) 



In a pure system Tjv(0)/2J = l/ln3. At small x 
T/v(x) ~ Tjv(0)(1 — A ! x) and, after some algebra, one 
obtains an analytical expression for A 1 



A 1 = - 



4 2 



3 In 3 



3 3 / 4 - 1 



- 1 



1.025 I 



+37 , (F10) 



which should be compared with the RPA answer A T RPA = 
4/3 |59| and an exact answer A I exact ~ +57 pq| ' One can 
see that in spite of the roughness of the approximation 
of the Ising spin degrees of freedom by bosons our ap- 
proach gives a good quantitative agreement with other 
approaches and an exact result. 
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